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Abstract 

We present a stochastic simulator for polycrystalline phase-change materials capable of spatio- 
temporal modelling of complex anneals. This is based on consideration of bulk and surface energies 
to generate rates of growth and decay of crystallites built up of 'monomers' that themselves may 
be quite complex molecules. We perform a number of simulations of this model using a Gillespie 
algorithm. The simulations are performed at molecular scale and using an approximation of local 
free energy changes that depend only on immediate neighbours. The sites are on a lattice chosen to 
have a lengthscale of the individual monomers, where each site gives information about a two-state 
local phase r (r = corresponds to amorphous and 1 corresponds to crystalline) and a continuous 
crystal orientation <p at each site. 

As an example we use this to model crystallisation in chalcogenide GST (GeSbTe) alloys used 
for example in phase-change memory devices, where reversible changes between amorphous and 
crystalline regimes are used to store and process information. We use our model to simulate anneals 
of GST including ones with non-trivial spatial and temporal variation of temperature; this gives 
good agreement to experimental incubation times at low temperatures while modelling non-trivial 
crystal size distributions and melting dynamics at higher temperatures. 

PACS numbers: 07.05.Tp (Computer modeling and Simulation), 64.60.De (Statistical mechanics of model 
systems; Ising model, Potts model, field-theory models, Monte Carlo techniques, etc). 
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I. INTRODUCTION 



This paper considers a model for phase change materials, i.e. alloys that can undergo 
reversible phase changes in response to anneals. Possible modelling techniques for poly- 
crystalline processes in reversible phase-change alloys such as GST (GeSbTe) range from 
molecular dynamic simulations at one end of the spectrum to empirical models at the other. 
The former are thermodynamically realistic but highly computer-intensive; the latter are 
fast, but hard to relate to material properties. Hence models used in practise tend to lie 



somewhere between the two, for example, Monte Carlo simulations [lOf] . JMAK and master 
equation based models 0, Q, Q or probabilistic cellular automata p, 0, Q • 

We describe the material as a 2D lattice of discrete 'sites' where each site is either 
crystalline or amorphous and there is an underlying orientation that varies continuously; 
these sites are on the lengthscale of monomers, though they do not correspond directly 
to individual monomers. As such our method for generating complex anneals, combines 
elements of probabilistic cellular automata (PCA) models P. fl?!. 18]. polycrystallinephase 
field models and uses thermodynamics from master equation models [9]. As 

we use a Gillespie algorithm for time-stepping we refer to it as a 'Phase Field Gillespie' 
model. In particular we retain a high level of simplicity because of discrete time and lattice 
space model, while retaining thermodynamic realism and hence keep fitting parameters to 
a minimum. 

Recall that crystallisation can be thought of as a two-stage process; nucleation (where a 
small crystallite needs to overcome an energy barrier dominated by interfacial energy) and 
growth (where the crystallite grows according to the availability of neighboring monomers 
and dominated by bulk energy). We assume each site has a set of locally-determined rate 
constants for transitions into a new state. These rates depend only on the current state of 
the site and that of its immediate neighbours. For the rates of growth and dissociation for 
GST we use thermodynamics from 13|, Il6j ]. 

The Gillespie algorithm 4j can be used to simulate the evolution under the assumption 
that the events are independent, instantaneous and never simultaneous. Each step of the 
algorithm has two parts; firstly it determines a random time to next event and secondly it 
determines which event occurs. This enables one to perform fast and physically plausible 
simulations of a number of crystallisation-related phenomena, including incomplete crystalli- 
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sation, melting and complex spatio-temporal anneals. As there are many possible events, 
one must efficiently use data-structures to ensure that the simulations run at a high speed 
and hence perform simulation of complex anneals in 2D on a standard desktop computer. 

After the algorithm for the Phase Field Gillespie simulator is described in Section [Til 
Section II II I presents the results of some bulk anneals of GST using this simulator, showing 
that the simulator can model nucleation effects, non-trivial anneals and melting. We include 
examples where the temperature depends on space and/or time; one can see a variety of 
effects and a good quantitative agreement with experimental temperature-dependent incu- 
bation times for GST. Finally Section IIVI discusses some possible extensions and limitations 
of the method. 



II. A PHASE FIELD GILLESPIE (PFG) CRYSTALLISATION SIMULATOR 

We consider a homogeneous (though not necessarily isotropic) material in 2D. The state 
of the material is described on a discrete regular lattice of grid points G on the length scale 
of the individual monomers. Each site is assumed to be either 'crystalline' or 'amorphous'. 
More precisely, at each grid point, (z, j) G G, the state is described by two quantities: 

• - a discrete 'phase' variable that is either (amorphous) or 1 (crystalline). 

• (pij - a continuous 'orientation' variable that varies over some range [0, ir) and gives a 
notional representation for the local orientation of the material. 

In particular we can determine that two adjacent crystalline sites are within the same crystal 
if and only if r i3 - = r H = 1 and 4>ij = 4>ki- 

The model we describe is a stochastic model that estimates rates of possible local changes 
to the state of the system (i.e. changes that affect only one site) and uses a Gillespie 



algorithm {4] to evolve the system in time. A Gillespie algorithm is optimal in that it will 
generate timesteps at a rate corresponding to the fastest rate that requires updating, though 
it is typically more complex to implement than Monte Carlo simulations [h]] . Although there 
are adaptations of the algorithm to other contexts [5j] we use the original version of Gillespie. 
We consider the following possible instantaneous events at a site (i, j) G G: 

• Nucleation - The site (i, j) and an adjacent site, originally both amorphous, become 



a crystallite at a rate Cf] 



nu 
ij ' 



• Growth - The site originally amorphous, becomes attached to an adjacent 
crystal of orientation ip at a rate CfL. 

• Dissociation - The site originally crystalline, detaches or dissociates from the 
crystal of which it is a part to become amorphous at a rate Cfj, and assumes a random 
orientation. 



A. The rate coefficients 

We approximate the rate coefficients for nucleation, growth and dissociation (C nu , C gr and 
C dl ) at each grid point by considering the change in bulk and surface energies of crystallites 
adjacent to that site. We define the set of neighbours of G G 

Nij = {{k,l) G G : (k, I) is a neighbour of and = \Nij\, 

the set of amorphous neighbours of 

N°™ = {(k, I) G N l3 : r kl = 0}, and <™ = |A* m |, 

and finally the set of crystalline neighbours of with a given orientation ijj 

AT- = {( k , I) g N t3 : cj> kl = $ and r« = 1}, and = \N^\. 

Note that <f, G {0, • • • , riyj. 

The rates are considered in a similar way to the derivation of master equation rates as 



in 13] and outlined below. We assume that 'interactions' between neighbours occur at a 
temperature-dependent rate 

R(T) = k Q et^) 

where E a is the activation energy and ks is the Boltzmann constant. The prefactor fco is 
used as a fitting parameter to normalise the results; see 
If adjacent sites have an 'interaction' we define 

rate at which a site transforms form amorphous to crystalline, 
resulting in a change A in surface of the crystallite. 

We assume local thermal equilibrium, meaning that the rate of the reverse transformation 
at an interaction is £ -1 (T, A). This rate varies with temperature as the bulk and surface 
energy vary. 



We compute the change in surface area of the crystallites on adding site to a neigh- 
bouring crystal of orientation ip by a linear approximation 



A — S n 



Tlij 



where S m is the surface area of a single site. This means that changing an isolated site in 
the middle of a crystal of orientation ip will result in a change A = —S m ( 
while creating a new crystal in the middle of a field of amorphous material will result in a 
change A = S m (as = 0). 

Putting this together (and noting that only by interaction with amorphous neighbours 
can a site nucleate) we get rate for nucleation that is 

( Eg \ ^am 

CT = i ko ^(T, S m ) if r v = 



if Tij = 1. 



The growth rate for an amorphous site to join a crystalline neighbour with orientation i[) is: 



CfL={ ^ ' n " ; 3 (2) 

1 if rij = 1. 

Finally, the dissociation rate for a crystalline site to become amorphous is: 

„,j I if r ii = 

1 fco eC"w) £(T, if r« = 1. 

B. The PFG Algorithm 

We now detail the operation of the PFG algorithm using the rates above. Initially, the 
whole domain is assumed to be an 'as deposited' amorphous state with a random distribution 
of 4>ij values and = 0, though one can restart the algorithm from any given state. 

For a square lattice we use eight neighbors, = 8, weighted according to their distance 
from the site. The new state of the site is then given by r'^ and using the stochastic 
simulation algorithm of Gillespie Q] as follows. This simulates up to a time T max . 



1. Start at time T = with given and 



2. Generate rate coefficients for all grid points C?" Cf-^ Cf- for nucleation, growth and 
dissociation respectively. We refer to these using a single index v = where 
a G {nu, (gr,V),di}. 



3. Compute the sum 



a 



^2 C iji> 



where ^ = {(fiti '■ (k,l) G iV^} is the set of orientations of neighbours to (i, j). 

4. Generate two independent random numbers r)i,r)2 uniformly distributed on (0, 1) and 
compute 

dr=-\og e (-). (4) 
Increment time to T = T + dr. If T > T max then stop. 

5. Identify the event v = a) corresponding to a G (nu, gr, di) and (k,l) with 
(pki = ip such that 

^2 a v < mao <^a u (5) 

6. Update the value of 0^- and r^. More precisely, perform the following updates accord- 
ing to corresponding reactions (nucleation, growth or dissociation) that occur: 

(a) Nucleation at (i, j); pick a (k,l) G N?™ 1 at random and set 

<P'ij = <t>'ki = 4>ij, ^ = r' kl = 1. 

(b) Growth from neighboring crystal with *0 = (fc, /) G iV^ into the amor- 
phous site set 

4>'ij = <t>kh r'^ = 1. 

(c) Dissociation at set 

rj,- = <^ = W, 

where is an independent random number uniformly distributed in the range 
of possible orientations [0, 7r). 

7. For the next iteration, copy 0^ = 4>'ij and = r^- and update the values of C° u Cf T 



ijip 



Of- 



8. Return to step 3 and recompute a . 
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Note that the main computational effort is the selection of the event (Step 5) based on 772; 
to minimize the number of operations needed to determine this we use a recursive bisection 
search and an efficient sorting of possible events. Also in the recomputation of rates (Step 
7) one can limit the updates to those sites that have changed and their neighbours. Finally, 
the computation of a (Step 3) in subsequent steps can considerably be accelerated by using 
only addition and subtraction of those rates that have changed. 



III. SIMULATIONS OF PHASE CHANGE FOR GST 



For the remainder of this paper we model the phase ch ange material GST used for 



read/write optical and electrical data storage devices, as in [13J]. Such a material has a 
fine balance between bulk and surface energies of crystals, meaning that one can find non- 
trivial nucleation and growth dynamics that varies with T. 

Let T m be the melting temperature; if we assume that the free energy change associated 
with crystallisation of a single site varies linearly with T — T m and the energy change asso- 
ciated with change in surface A is a A with a constant, then the rate £(T, A) can be written 
as 

£(T,A) =exp 



L I 1.0 - X-^ ( ° A 



Following [13j], we assume that 



(6) 



(7) 



where constants are a = 2.2 x 10~ 6 Jcm~ 2 , the interfacial energy density between amor- 
phous and crystalline phases S m = 2.1187 x 10~ 14 cm 2 the molecular surface area of the 
material. We use E a = 2.1 eV and k® = 10 16 /is _1 where we use time units of microseconds 
for convenience. The other constants are as follows: 

• AHf = 625 Jcm~ 3 is the enthalpy of fusion from the data obtained from differential 
scanning calorimeter experiments on GST. 

• v m — 2.9 x 10~ 22 cm 3 is the molecular volume of GST and S m = 2.1 x 10~ 14 cm 2 
assuming approximately spherical shape. 

• 7^=889°^ is the melting temperature 

• k B = 1.381 x 10~ 23 J/°K is Boltzmann's constant. 



Using these values we obtain L = 7.3816. In our simulations we assume we have N 2 sites 
with periodic boundary conditions applied in both directions; i.e. r i+N j = Tij +N = r^. The 
parameters for the Phase Field Gillespie algorithm outlined above give realistic quantitative 
agreement with crystal growth in GST over a range of temperatures. 



A. Nucleation and crystal growth 

We simulate using an N 2 grid with N = 256. Note that the crystalline fraction X for 
such a grid can be calculated as 

x = E r v 

(ij)eG 

where clearly < X < 1 and X = 1 corresponds to a fully crystalline state. 

We show in Figure [1] the increase in the crystalline fraction X as a function of time 
starting at fully amorphous for T = 131°C; after an initial incubation the fraction quickly 
increases to saturate near X = 1. The insets show that the growth occurs subject to 
random fluctuations from the algorithm. Near X = 1 there is still a nontrivial process of 
detachment and reattachment of sites from crystals that leads to grain coarsening over a long 
timescale. Figure [2] shows the progress of this anneal at three stages; soon after inception, 
at approximately 20% progress and in a polycrystalline state, while Figure [3] shows the 
development of the distribution of crystal sizes as the anneal progresses. 

The incubation time (defined here as the time to 20% crystallinity from fully amorphous) 
is shown against temperature in Figure H] and for comparison the results from data from 
experiments [jjj] as well as for the master equation model [ij are plotted. We note that the 
Phase Field Gillespie simulations show a temperature dependence that is close to experi- 
mental results of 15J both in form and value. As with the master equation model, there 
is effectively only one fitting parameter in the model, the prefactor k Q and this is fixed 
independent of temperature. 



B. Spatio-temporal anneals 

One can easily adapt the algorithm to the case where the temperature, and therefore the 
rates of the reactions, depends on the spatial location; the algorithm is exactly as presented 
before except that T now depends Tij on site and time. As an example, in Figure [5] shows 
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the development of a band of GST material that is held at 227°C on the left boundary and 
A77°C on the right boundary for a period of time. On the left hand side the growth is very 
slow while on the right the nucleation energy is more difficult to overcome, meaning that 
initial growth is fastest in the intermediate region. A final example is given in Figure [6] 
where a sample is subjected to a complex sequence of spatio-temporal anneals; see caption 
for details. 



IV. DISCUSSION 



The Phase Field Gillespie algorithm introduced in the paper incorporates features from 
a few different models of crystallisation, and can be thought of as a thermodynamically 
motivated caricature of a molecular simulation. We highlight a few features that could be 
investigated to make more physically realistic models: 

• The current model is based on a 2D grid meaning that the interpretation of the volume 
and surface area of the monomers at each site should depend on interfacial energies, 
or alternatively this could be adapted to a 3D grid with suitable boundary conditions. 

• We assume the energies of the crystallites do not depend on orientation. It would be 
relatively easy to include anisotropy, meaning that crystallites should grow depending 
on their orientation. 

• We have so far only considered the behavior of the crystallisation by imposing a 
temperature that may be uniform, or may have spatial non-uniformities. It would 
be of interest to investigate the coupling of this to phase, for example as occurs in 
electrical heating of GST. In this case onset of percolation results in lower resistance 
and hence increase heating. 

Nonetheless the current model can evidently produce reasonably realistic and numerically ef- 
ficient simulations of crystallisation behaviour even for complex spatio-temporal anneals and 
as such we believe the model is worthy of further investigation. We also suggest that these 
techniques will be useful for modelling phase change devices that use reversible transitions 



in GST alloys to perform computations and for multi-state storage 
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FIG. 1: Crystalline fraction X as a function of time during low temperature anneals at 131°C. 
Detail of the progress of the anneal is shown during the growth phase and when the crystalline 
fraction has saturated near X = 1. 
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FIG. 2: Images showing progress in crystallisation (a-c) for T = 131°C (as shown in Figure [T]); and 
(d) for T = 407°C, starting with pure amorphous material. The colours are assigned arbitrarily 
to different oriented crystal grains, (a) shows after 2000 steps of the algorithm a number of nuclei 
with X = 0.0507 after time 248s, (b) shows after 10000 steps with X = 0.204 after time 743s, (c) 
shows after 10 5 steps with X = 0.999115 after time 68930s. Similarly, (d) shows the state after 10 6 
steps corresponding to 4.347/zs at the higher temperature. Observe the faster progress and larger 
crystals that result at the higher temperature. 
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FIG. 3: The relative frequency of crystallites of different sizes corresponding to the (a), (b) and 
(c) of Figure El Observe the peak in crystal size distribution at size 15-20 for the developed crystal 
structure (c). 
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FIG. 4: Incubation times given by Phase Field Gillespie simulations, with master equation simula- 
tions of GST crystallisation from [1] and experimental data from [15J shown for comparison. Note 
that the Phase Field Gillespie simulation produces an excellent agreement with experiment. 
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FIG. 5: Images showing progress in crystallisation for a sample held in a temperature gradient 
where the left boundary is 227°C and the right is 477°C. Observe the appearence of a band of 
higher crystallinity as time progresses from (a) after 17.6ns, (b) after 70ns, (c) after 554ns and 
(d) after 22.9//S. Observe that the effective nucleation size is larger on the right (hotter) side of 
the sample. 
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FIG. 6: Progress of a multi-step anneal, demonstrating both spatial and temporal variation in tem- 
peratures. Starting from amorphous, a sample is first subjected lus of a linear spatial temperature 
gradient, the left at 227° C and the right at 477°C; (a) shows resulting the crystal structure. For 
the next 0.1s it is maintained at 227°C and in doing so progresses towards almost complete crys- 
tallisation but with a clear banded structure shown in (b). Finally the sample is raised to 477° C 
for 15ns which is enough for the crystals to almost entirely dissociate; the structure at X = 0.5 is 
shown in (c) and the final state is (d). Although below melting temperature, the critical nucleus 
size is too large for crystals of this size to survive, (e) Shows the crystalline fraction X as a function 
of algorithm step; the timesteps vary by many orders of magnitude as the anneal progresses. 
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